disp('tjalff');

T = 0.1;
H = 0.0;
h0 = H/T;
mumu = 1.0;
Nc = 4;

Num = 400;

vk = linspace(0, 4, Num);
vep_ori = vk.^2 - mumu - T*log(2*cosh(h0/2));


F = @(x, T) -T*log(1+exp(-x/T));
F2 = @(u, m) log((1 + exp(-u)*fn_inveta(m, h0))...
    ./(1 + fn_inveta(m, h0)));

% please notice that the mu is double used
mu = zeros(Nc, length(vk));

vep = vep_ori;

while 1
    
    vep0 = vep;
    mu0 = mu;
    %%%%%%%% for epsilon  %%%%%%%%%%%%%%%%%%%%%%%%%%
    vep = vep_ori + fn_anf(2, vk, F(vep,T), vk);
    for n = 1:Nc
        vu = mu(n,:);
        vep = vep - T*fn_anf(n, vk, F2(vu,n), vk);
    end
    %%%%%%%% for string  %%%%%%%%%%%%%%%%%%%%%%%%%%
    for n = 1:Nc
        vu = -1/T*fn_anf(n, vk, F(vep, T), vk);
        for m = 1:Nc
            vu = vu + fn_Tnm(n,m,fn_inveta(m, h0),vk,mu(m,:),vk);
        end
        mu(n,:) = vu;
    end
    
    fprintf('diss = %e   ', fn_norm(vep - vep0));
    fprintf('diss_m = %e \n', fn_norm(mu - mu0));
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if fn_norm(vep - vep0)<1e-8 && fn_norm(mu - mu0)<1e-8
        break;
    end
    
end

close(figure(4));
figure(4); hold on
for n = 1:Nc
    plot(fn_ex1(vk), fn_ex2(mu(n,:)), 'r-*');
end

vx_nt = fn_nt(vk);
pressure = -2*sum(vx_nt.*F(vep,T))/2/pi;
fprintf('the pressure = %8.8f \n', pressure);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Num_point = 400;

vNc = 1:4;
vpressure = zeros(1, length(vNc));
for i = 1:length(vNc)
    Nc = vNc(i);
    [vk, vep, ck, cX] = fn_solve_tba(T,H,mumu,Nc,Num_point);
    vx_nt = fn_nt(vk);
    pressure = -2*sum(vx_nt.*F(vep,T))/2/pi;
    fprintf('Nc = %d, the pressure = %8.8f \n', Nc, pressure);
end

close(figure(9));
figure(9);
plot(vNc, vpressure, '-ro');


figure(24);
plot(fn_ex1(vk), fn_ex2(vep), '-bo');

figure(4); 
for n = 1:Nc
    
    figure(4); hold on;
    plot(fn_ex1(ck{n}), fn_ex2(fn_X2u(n, h0, cX{n})), 'bo');
    xlim([-10,10]);
    
end

figure(14); 
for n = 1:Nc
    
    figure(14); hold on;
    plot(fn_ex1(ck{n}), fn_ex2(cX{n}), '-bo');
    xlim([-20,20]);
    
end


disp('tjsjlftj');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [vk, vep, ck, cX] = fn_solve_tba(T,H,mumu,Nc,Num_point)

h0 = H/T;
Num = Num_point;

vk = linspace(0, 3*sqrt(mumu + T*log(2*cosh(h0/2))), Num);
vep_ori = vk.^2 - mumu - T*log(2*cosh(h0/2));

F = @(x, T) -T*log(1+exp(-x/T));
% F2 = @(u, m) log((1 + exp(-u)*fn_inveta(m, h0))...
%     ./(1 + fn_inveta(m, h0)));
ck = cell(1, Nc);
cX = cell(1, Nc);

for n = 1:Nc
    ck{n} = linspace(0, 4, Num);
    cX{n} = zeros(1, Num);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vep = vep_ori;
while 1
    
    vep0 = vep;
    ck0 = ck;
    cX0 = cX;
    
    vep = vep_ori + fn_anf(2, vk, F(vep,T), vk);
    
    for n = 1:Nc
        vep = vep + T*fn_anf(n, ck{n}, cX{n}, vk);
    end
    
    vnorm = zeros(1, Nc);
    for n = 1:Nc
        
        vFT = -1/T*F(vep, T);
        [vp, vr] = fn_iteration(n,h0,vk, vFT, ck, cX, Num);
        ck{n} = vp;
        cX{n} = vr;
        
        vnorm(n) = fn_norm_2(ck{n},cX{n},ck0{n},cX0{n});
        
    end
    
    fprintf('diss = %1.2e   ', fn_norm(vep - vep0));
    fprintf('diss = %1.2e   \n', fn_norm(vnorm)); 
    
     if fn_norm(vep - vep0)<1e-8 && fn_norm(vnorm)<1e-8
        break;
    end

end

end


function res = fn_norm_2(vx1, vf1, vx2, vf2)

vx = [vx1, vx2];

x1 = vx1(end);
x2 = vx2(end);

xmin = min(x1, x2);
vx(vx>xmin) = [];
vx = sort(vx);

vf1_ins = interp1(vx1, vf1, vx);
vf2_ins = interp1(vx2, vf2, vx);

res = fn_norm(vf1_ins-vf2_ins);

end



function vX = fn_u2X(m, h0, vu)

a = fn_inveta(m, h0);
vX = -log((1 + a*exp(-vu))./(1 + a));

end


function vu = fn_X2u(m, h0, vX)

a = fn_inveta(m, h0);
vu = -log(1/a*(exp(-vX)-1) + exp(-vX));

end


function res = fn_f(n, k, a, vk, vFT, ck, cX)

Num = length(cX);
res = fn_int_an(n, vk, vFT, k);

for m = 1:Num
    
   res = res - fn_int_Tnm(n, m, ck{m}, cX{m}, k); 
    
end

res = -log((1 + exp(-res)*a)/(1+a));

end


function [vp, vr] = fn_iteration(n,h0,vk,vFT,ck,cX,Num_point)

a = fn_inveta(n,h0);
f = @(x) fn_f(n, x, a, vk, vFT, ck, cX); 

p_left = 0;
p_right = 1000;

vp = [p_left, p_right];
vr = [f(p_left), f(p_right)];

while 1
    
    vd = abs(vr(2:end) - vr(1:end-1));
    
    [~, max_region] = max(vd);
    
    ind_left = max_region;
    ind_right = max_region + 1;
    
    p_left = vp(ind_left);
    p_right = vp(ind_right);
    p_mid = 1/2*(p_left + p_right);
    
    new_r = f(p_mid);
    
    vp = [vp(1:ind_left), p_mid, vp(ind_right:end)];
    vr = [vr(1:ind_left), new_r, vr(ind_right:end)];

    if length(vp) >= Num_point
        break;
    end

end

end



function [vp, vr] = fn_generate_curve(n, Num_point, vx, vf)

f = @(x) fn_int_an(n, vx, vf, x);

p_left = 0;
p_right = 10000;

vp = [p_left, p_right];
vr = [f(p_left), f(p_right)];

while 1
    
    vd = abs(vr(2:end) - vr(1:end-1));
    
    [~, max_region] = max(vd);
    
    ind_left = max_region;
    ind_right = max_region + 1;
    
    p_left = vp(ind_left);
    p_right = vp(ind_right);
    p_mid = 1/2*(p_left + p_right);
    
    new_r = f(p_mid);
    
    vp = [vp(1:ind_left), p_mid, vp(ind_right:end)];
    vr = [vr(1:ind_left), new_r, vr(ind_right:end)];

    if length(vp) >= Num_point
        break;
    end

end

end


function res = fn_int_Tnm(n, m, vx, vf, k)

res = 0;

if n == m
    
    for i = 2:2:(2*n-2)
        
       res = res + 2*fn_int_an(i, vx, vf, k);  
        
    end
    
    res = res + fn_int_an(2*n, vx, vf, k);
    
else
    
    for i = (abs(m-n)+2):2:(m+n-2)
        
        res = res + 2*fn_int_an(i, vx, vf, k);
        
    end
    
    res = res + fn_int_an(abs(m-n), vx, vf, k);
    res = res + fn_int_an(m+n, vx, vf, k);
    
end

end


function res = fn_int_an(n, vx, vf, k)

res = fn_int_a1(vx/n, vf, k/n);

end


function res = fn_int_a1(vx, vf, k)

vp_tail = vx(end-9:end);
vf_tail = vf(end-9:end);
v2 = 1./vp_tail.^2;
v4 = 1./vp_tail.^4;
kc = vx(end);
[alpha, beta] = fn_fit(v2, v4, vf_tail);

vp0 = linspace(-4, 4, 400);
vp0 = vp0.^2.*sign(vp0);
vp1 = vp0 + k;
vp2 = vx;
vp3 = vp0 - k;

vp = [vp1, vp2, vp3];
vp(vp<0) = [];
vp(vp>kc) = [];
vp = sort(vp);

vf2 = interp1(vx, vf, vp);
va = 1/2/pi./(1/4 + (k - vp).^2) + 1/2/pi./(1/4 + (k + vp).^2);

vnt = fn_nt(vp);
res = sum(vnt.*va.*vf2);

res = res + alpha*fn_ff2(k, kc) + beta*fn_ff4(k, kc);

end



function [v_alpha, vy_fit] = fn_fit_an(vx, vy, Num)

an = @(x,n) 1/2/pi*n./((n/2).^2 + x.^2);

vx = transpose(vx);
vy = transpose(vy);
vy_fit = 0;

c_an = cell(1, Num);
for n = 1:Num
    
   c_an{n} = an(vx, n); 
    
end

V = zeros(Num, 1);
for n = 1:Num
    
   V(n) = transpose(c_an{n})*vy; 
    
end

M = zeros(Num, Num);
for q = 1:Num
    for p = 1:Num
        
        M(q, p) = transpose(c_an{q})*c_an{p};
        
    end
end

v_alpha = transpose(M\V);
fprintf('condition = %f \n', cond(M));

for n = 1:Num
    
   vy_fit = vy_fit + v_alpha(n)*c_an{n}; 
    
end

vy_fit = transpose(vy_fit);

end



function res = fn_test(h)

vn = 1:1000;
res = sum(log(1 - sinh(h/2).^2./sinh(h/2*(vn+1)).^2));

end



function res = fn_G(x, T)

if x/T>33
    
    res = 0;
    
else

    res = 1./(1 + exp(x/T));
    
end

end



function Pn = fn_Pn(H, T, un)

[Nc, Num] = size(un);
Pn = zeros(Nc, Num);

if H>0

    for n = 1:Nc
        
        vu = un(n,:);
        P = (1 - exp(-vu))/T.*((n+1)*coth(H/2/T*(n+1)) - coth(H/2/T)) ...
            .*2.*sinh(H/2/T)^2 ...
            ./(cosh(H/2/T*(2*n+2)) - exp(-vu) + (exp(-vu) - 1).*cosh(H/2/T*2));
        Pn(n,:) = P;
        
    end

end

end

function Qn = fn_Qn(H, T, un)

[Nc, Num] = size(un);
Qn = zeros(Nc, Num);
h0 = H/T;

for n = 1:Nc
    vu = un(n,:);
    eta = fn_eta(n, h0);  
    Q = 1./(eta*exp(vu) + 1);
    Qn(n,:) = Q;
end

end



function [pressure_dh, vx, vxn, vep_dh, mu_dh] = fn_pressure_dh(mumu, H, T, Nc)

h0 = H/T;
F = @(x, T) -T*log(1+exp(-x/T));
F2 = @(u, m) log((1 + exp(-u)*fn_inveta(m, h0))...
    ./(1 + fn_inveta(m, h0)));

Num = 400;
ep_cut = 2*sqrt(mumu);
u_cut = (1:Nc)*10;

vx = linspace(0, ep_cut, Num);
vxn = zeros(Nc, Num);

for n = 1:Nc
    vxn(n,:) = linspace(0,u_cut(n),Num);
end

vep = vx.^2 - mumu - T*log(2*cosh(H/2/T));
vep_ori = vep;
% please notice that the mu is double used
mu = zeros(Nc, Num);
while 1
    
    vep0 = vep;
    mu0 = mu;
    
    %%%%%%%% for epsilon  %%%%%%%%%%%%%%%%%%%%%%%%%%
    vep = vep_ori + fn_anf(2, vx, F(vep,T), vx);   
    for n = 1:Nc     
        vu = mu(n,:);
        vep = vep - T*fn_anf(n, vxn(n,:), F2(vu,n), vx);     
    end
    
    %%%%%%%% for string  %%%%%%%%%%%%%%%%%%%%%%%%%%
    for n = 1:Nc
        vu = -T*fn_anf(n, vx, F(vep, T), vxn(n,:));
        for m = 1:Nc       
            vu = vu + fn_Tnm(n,m,fn_inveta(m, h0),vxn(m,:),mu(m,:),vxn(n,:));       
        end
        mu(n,:) = vu;
    end
    
    fprintf('diss = %e   ', fn_norm(vep - vep0));
    fprintf('diss_m = %e \n', fn_norm(mu - mu0));
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    if fn_norm(vep - vep0)<1e-8 && fn_norm(mu - mu0)<1e-8
        break;
    end
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x1 = fn_find_zero(vx, vep + 16*T);
x2 = fn_find_zero(vx, vep + 8*T);
% x0 = fn_find_zero(vx, vep);
x3 = fn_find_zero(vx, vep - 8*T);
x4 = fn_find_zero(vx, vep - 16*T);


while 1
    
    x_left = x2;
    x_right = x3;
    
    Num2 = 200;
    vx_left = linspace(0, x1, Num2);
    vx_mid = linspace(x1, x4 ,Num2);
    vx_right = linspace(x4, ep_cut, Num2);
    vx = [vx_left,vx_mid(2:end),vx_right(2:end)];
    
    vxn = zeros(Nc, Num);
    for n = 1:Nc
        vxn(n,:) = linspace(0,u_cut(n),Num);
    end
    
    vep = vx.^2 - mumu - T*log(2*cosh(H/2/T));
    vep_ori = vep;
    % please notice that the mu is double used
    mu = zeros(Nc, Num);
    while 1
        
        vep0 = vep;
        mu0 = mu;
        
        %%%%%%%% for epsilon  %%%%%%%%%%%%%%%%%%%%%%%%%%
        vep = vep_ori + fn_anf(2, vx, F(vep,T), vx);
        for n = 1:Nc
            vu = mu(n,:);
            vep = vep - T*fn_anf(n, vxn(n,:), F2(vu,n), vx);
        end
        
        %%%%%%%% for string  %%%%%%%%%%%%%%%%%%%%%%%%%%
        for n = 1:Nc
            vu = -T*fn_anf(n, vx, F(vep, T), vxn(n,:));
            for m = 1:Nc
                vu = vu + fn_Tnm(n,m,fn_inveta(m, h0),vxn(m,:),mu(m,:),vxn(n,:));
            end
            mu(n,:) = vu;
        end
        
        fprintf('diss = %e   ', fn_norm(vep - vep0));
        fprintf('diss_m = %e \n', fn_norm(mu - mu0));
        fprintf(' Nc = %d \n', Nc);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if fn_norm(vep - vep0)<1e-8 && fn_norm(mu - mu0)<1e-8
            break;
        end
        
    end
    
    x1 = fn_find_zero(vx, vep + 16*T);
    x2 = fn_find_zero(vx, vep + 8*T);
    x0 = fn_find_zero(vx, vep);
    x3 = fn_find_zero(vx, vep - 8*T);
    x4 = fn_find_zero(vx, vep - 16*T);
    
    if x_left<x0 && x_right>x0
        break;
    end
    
end

Pn = fn_Pn(H, T, mu);
Qn = fn_Qn(H, T, mu);

vep_dh_ori = -1/2*tanh(H/2/T)*ones(1, length(vx));
vep_dh = vep_dh_ori;
mu_dh = zeros(Nc, Num);

while 1
    
    vep_dh0 = vep_dh;
    mu_dh0 = mu_dh;
    
    %%%%%%%% for epsilon  %%%%%%%%%%%%%%%%%%%%%%%%%%
    vep_dh = vep_dh_ori + fn_anf(2, vx, fn_G(vep,T).*vep_dh, vx);   
    for n = 1:Nc     
        vep_dh = vep_dh - T*fn_anf(n, vxn(n,:),Pn(n,:)-Qn(n,:).*mu_dh(n,:), vx);     
    end
    
    %%%%%%%% for string  %%%%%%%%%%%%%%%%%%%%%%%%%%
    for n = 1:Nc
        vu = -T*fn_anf(n, vx, fn_G(vep,T).*vep_dh, vxn(n,:));
        for m = 1:Nc       
            vu = vu + fn_Tnm0(n,m,vxn(m,:),Pn(m,:)-Qn(m,:).*mu_dh(m,:),vxn(n,:));       
        end
        mu_dh(n,:) = vu;
    end
    
    fprintf('diss_dh = %e   ', fn_norm(vep_dh - vep_dh0));
    fprintf('diss_m_dh = %e \n', fn_norm(mu_dh - mu_dh0));
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    if fn_norm(vep_dh - vep_dh0)<1e-8 && fn_norm(mu_dh - mu_dh0)<1e-8
        break;
    end
    
end

vx_nt = fn_nt(vx);
pressure_dh = -2*sum(vx_nt.*vep_dh.*fn_G(vep, T))/2/pi;

end



function [pressure, vx, vxn, vep, mu] = fn_pressure(mumu, H, T, Nc)

h0 = H/T;
F = @(x, T) -T*log(1+exp(-x/T));
F2 = @(u, m) log((1 + exp(-u)*fn_inveta(m, h0))...
    ./(1 + fn_inveta(m, h0)));

Num = 400;
ep_cut = 2*sqrt(mumu);
u_cut = (1:Nc)*10;

vx = linspace(0, ep_cut, Num);
vxn = zeros(Nc, Num);

for n = 1:Nc
    vxn(n,:) = linspace(0,u_cut(n),Num);
end

vep = vx.^2 - mumu - T*log(2*cosh(H/2/T));
vep_ori = vep;
% please notice that the mu is double used
mu = zeros(Nc, Num);
while 1
    
    vep0 = vep;
    mu0 = mu;
    
    %%%%%%%% for epsilon  %%%%%%%%%%%%%%%%%%%%%%%%%%
    vep = vep_ori + fn_anf(2, vx, F(vep,T), vx);   
    for n = 1:Nc     
        vu = mu(n,:);
        vep = vep - T*fn_anf(n, vxn(n,:), F2(vu,n), vx);     
    end
    
    %%%%%%%% for string  %%%%%%%%%%%%%%%%%%%%%%%%%%
    for n = 1:Nc
        
        % changed !!!!!!!!!!!!!!!!!!!!!!
        vu = -1/T*fn_anf(n, vx, F(vep, T), vxn(n,:));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        for m = 1:Nc       
            vu = vu + fn_Tnm(n,m,fn_inveta(m, h0),vxn(m,:),mu(m,:),vxn(n,:));       
        end
        mu(n,:) = vu;
    end
    
    fprintf('diss = %e   ', fn_norm(vep - vep0));
    fprintf('diss_m = %e \n', fn_norm(mu - mu0));
    fprintf(' Nc = %d \n', Nc);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    if fn_norm(vep - vep0)<1e-8 && fn_norm(mu - mu0)<1e-8
        break;
    end
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x1 = fn_find_zero(vx, vep + 16*T);
x2 = fn_find_zero(vx, vep + 8*T);
x0 = fn_find_zero(vx, vep);
x3 = fn_find_zero(vx, vep - 8*T);
x4 = fn_find_zero(vx, vep - 16*T);

% close(figure(9));
% figure(9); hold on;
% plot(vx, vep, '-bo');
% plot(x0, 0, 'g+');

while 1
    
    x_left = x2;
    x_right = x3;
    
    Num2 = 200;
    vx_left = linspace(0, x1, Num2);
    vx_mid = linspace(x1, x4 ,Num2);
    vx_right = linspace(x4, ep_cut, Num2);
    vx = [vx_left,vx_mid(2:end),vx_right(2:end)];
    
    vxn = zeros(Nc, Num);
    for n = 1:Nc
        vxn(n,:) = linspace(0,u_cut(n),Num);
    end
    
    vep = vx.^2 - mumu - T*log(2*cosh(H/2/T));
    vep_ori = vep;
    % please notice that the mu is double used
    mu = zeros(Nc, Num);
    while 1
        
        vep0 = vep;
        mu0 = mu;
        
        %%%%%%%% for epsilon  %%%%%%%%%%%%%%%%%%%%%%%%%%
        vep = vep_ori + fn_anf(2, vx, F(vep,T), vx);
        for n = 1:Nc
            vu = mu(n,:);
            vep = vep - T*fn_anf(n, vxn(n,:), F2(vu,n), vx);
        end
        
        %%%%%%%% for string  %%%%%%%%%%%%%%%%%%%%%%%%%%
        for n = 1:Nc
            
            % changed !!!!!!!!!!!!!!!!!!!!!!
            vu = -1/T*fn_anf(n, vx, F(vep, T), vxn(n,:));
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
            for m = 1:Nc
                vu = vu + fn_Tnm(n,m,fn_inveta(m, h0),vxn(m,:),mu(m,:),vxn(n,:));
            end
            mu(n,:) = vu;
        end
        
        fprintf('diss = %e   ', fn_norm(vep - vep0));
        fprintf('diss_m = %e \n', fn_norm(mu - mu0));
        fprintf(' Nc = %d \n', Nc);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if fn_norm(vep - vep0)<1e-8 && fn_norm(mu - mu0)<1e-8
            break;
        end
        
    end
    
    
    x1 = fn_find_zero(vx, vep + 16*T);
    x2 = fn_find_zero(vx, vep + 8*T);
    x0 = fn_find_zero(vx, vep);
    x3 = fn_find_zero(vx, vep - 8*T);
    x4 = fn_find_zero(vx, vep - 16*T);
    
    if x_left<x0 && x_right>x0
        break;
    end
    
end

% figure(9); hold on;
% plot(vx, vep, '-ro');
% plot(x0, 0, 'gs');

vx_nt = fn_nt(vx);
pressure = -2*sum(vx_nt.*F(vep,T))/2/pi;

end


function res = fn_int_with_tail(vx, vf)

vx_nt = zeros(1,length(vx));
vx_len = vx(2:end) - vx(1:end-1);
vx_nt(1:end-1) = 1/2*vx_len;
vx_nt(2:end) = vx_nt(2:end) + 1/2*vx_len;

res =sum(vx_nt.*vf);

kc = vx(end);
v_tail = vx(end-9:end);
v2 = 1./v_tail.^2;
v4 = 1./v_tail.^4;
vf_tail = vf(end-9:end);
[alpha, beta] = fn_fit(v2, v4, vf_tail);

res = res + alpha/kc + beta/3/kc^3;

res = res*2;

end



function res = fn_yn_inf(n, h0)

if h0>0
    
    res = n*h0 ...
        + log((1+exp(-(2*n+2)*h0)-exp(-n*h0)-exp(-(n+2)*h0)) ...
        ./(1+exp(-2*h0)-2*exp(-h0)));

else
   
    res = log(n.^2 + 2*n);
    
end

end


function vu_out = fn_Tnm(n, m, a, vx, vu, vx_out)

vu_out = 0;

vf = log((1 + exp(-vu)*a)./(1 + a));

if m == n
    
    for i = 2:2:(2*n-2)
        
        vf_out = fn_anf(i, vx, vf, vx_out);
        vu_out = vu_out + 2*vf_out;
        
    end
    
    vu_out = vu_out + fn_anf(2*n, vx, vf, vx_out);
    
else
    
    vu_out = vu_out + fn_anf(abs(m-n), vx, vf, vx_out);
    vu_out = vu_out + fn_anf(m+n, vx, vf, vx_out);
    
    for i = (abs(m-n)+2):2:(m+n-2)
        
        vu_out = vu_out + 2*fn_anf(i, vx, vf, vx_out);
        
    end
    
end

end



function vu_out = fn_Tnm0(n, m, vx, vf, vx_out)

vu_out = 0;

if m == n
    
    for i = 2:2:(2*n-2)
        
        vf_out = fn_anf(i, vx, vf, vx_out);
        vu_out = vu_out + 2*vf_out;
        
    end
    
    vu_out = vu_out + fn_anf(2*n, vx, vf, vx_out);
    
else
    
    vu_out = vu_out + fn_anf(abs(m-n), vx, vf, vx_out);
    vu_out = vu_out + fn_anf(m+n, vx, vf, vx_out);
    
    for i = (abs(m-n)+2):2:(m+n-2)
        
        vu_out = vu_out + 2*fn_anf(i, vx, vf, vx_out);
        
    end
    
end

end




function res = fn_inveta(m, h0)

if h0>0
    res = sinh(h0/2)^2/(sinh(h0/2*(m+1))^2 - sinh(h0/2)^2);
else
    res = 1/(m^2 + 2*m);
end
     
end


function res = fn_eta(m, h0)

if h0>0
    res = (sinh(h0/2*(m+1))^2 - sinh(h0/2)^2)/sinh(h0/2)^2;
else
    res = m^2 + 2*m;
end

end



function vf_out4 = fn_anf(n, vx, vf, vx_out)

an = @(x) 1/2/pi*n./((n/2).^2 + x.^2);

vk = vx_out;
vp = vx;

mk = transpose(vk)*ones(1,length(vp));
mp = ones(length(vk),1)*vp;
ma = an(mk - mp) + an(mk + mp);

vpd = vp(2:end) - vp(1:end-1);
vp_nt = zeros(1,length(vp));
vp_nt(1:end-1) = 1/2*vpd;
vp_nt(2:end) = vp_nt(2:end) + 1/2*vpd;

vf_out = ma*diag(vp_nt)*transpose(vf);
vf_out = transpose(vf_out);
vf_out0 = vf_out;

kc = vp(end);
v_tail = vp(end-9:end);
v2 = 1./v_tail.^2;
v4 = 1./v_tail.^4;
vf_tail = vf(end-9:end);
[alpha, beta] = fn_fit(v2, v4, vf_tail);

%  fprintf('alpha = %f, beta = %f \n', alpha, beta);

vf_out_addon2 = alpha ...
    *8/2/pi/n^2*(fn_f2(2*vk/n, 2*kc/n) + fn_f2(-2*vk/n, 2*kc/n));

vf_out_addon4 = beta ...
    *32/2/pi/n^4*(fn_f4(2*vk/n, 2*kc/n) + fn_f4(-2*vk/n, 2*kc/n));

%  vf_out2 = vf_out0 + vf_out_addon2;

vf_out4 = vf_out0 + vf_out_addon2 + vf_out_addon4;

end



function [alpha, beta] = fn_alpha_and_beta(vx, vf)

vx_tail = vx(end-9:end);
vf_tail = vf(end-9:end);
v2 = 1./vx_tail.^2;
v4 = 1./vx_tail.^4;

[alpha, beta] = fn_fit(v2, v4, vf_tail);

end


function vx_ex1 = fn_ex1(vx)

vx_ex1 = [-vx(end:-1:2), vx];

end

function vx_ex2 = fn_ex2(vx)

vx_ex2 = [vx(end:-1:2), vx];

end


function res = fn_int_tail(k, kc, n)

res = 8/pi*n*log(4)*k./(4*k.^2 + n^2).^2 ...
    + (4*k.^2 - n^2)*2./(4*k.^2 + n^2).^2 ...
    + 2/pi*n./(4*k.^2 + n^2)./kc ...
    - 4/pi*(4*k.^2 - n^2)./(4*k.^2 + n^2).^2.*atan(2*(k + kc)./n) ...
    + 8/pi*n*k./(4*k.^2 + n^2).^2.*log(kc.^2./(n^2 + 4*(k + kc).^2));

end


function res = fn_ff2(k, kc)

res = 2/pi/kc/(4*k^2+1) ...
    +4/pi*(4*k^2-1)/(4*k^2+1)^2*(pi/2+atan(2*(k-kc))) ...
    +8*k/pi/(4*k^2+1)^2*log((1+4*(k-kc)^2)/(4*kc^2));

end

function res = fn_ff4(k, kc)

res = 2/3/pi/(1+4*k^2)/kc^3 ...
    +8*k/pi/(1+4*k^2)^2/kc^2 ...
    +8*(48*k^4+8*k^2-1)/pi/(1+4*k^2)^4/kc ...
    +16/pi*(16*k^4-24*k^2+1)/(1+4*k^2)^4*(pi/2+atan(2*(k-kc))) ...
    +64*k/pi*(4*k^2-1)/(4*k^2+1)^4*log((1+4*(k-kc)^2)/(4*kc^2));

end


function res = fn_f2(x, kc)

res = -pi/2*(1 - x.^2)./(1 + x.^2).^2 ...
    + 1/kc./(1 + x.^2) + (1 - x.^2)./(1 + x.^2).^2.*atan(kc - x) ...
    + x./(1 + x.^2).^2.*log((1 + (kc - x).^2)./kc^2);

end


function res = fn_f4(x, kc)

res = 3*pi*(1 - 6*x.^2 + x.^4) ...
    + 6*(1 - 6*x.^2 + x.^4).*atan(x - kc) ...
    + 12*x.*(x.^2 - 1).*log((1 + (x - kc).^2)./kc^2) ...
    + 2*(1 + x.^2).^3/kc^3 ...
    + 6*x.*(1 + x.^2).^2/kc^2 ...
    + 6*(1 + x.^2).*(3*x.^2 - 1)/kc;

res = res./6./(1 + x.^2).^4;

end



function x0 = fn_find_zero(vx, vep)

[a0, a2, a4] = fn_fit_poly(vx, vep);
x0 = sqrt((sqrt(a2^2-4*a0*a4)-a2)/2/a4);

end



function [a0, a2, a4] = fn_fit_poly(vx, vep)

n_right = find(vep>0, 1);
n_left = n_right - 1;

vind = [(n_left - 4):n_left, n_right:(n_right + 4)];

x = vx(vind);
y = vep(vind);

Num = length(x);

V0 = ones(1,Num);
V2 = x.^2;
V4 = x.^4;

n = @(Vec1, Vec2) Vec1*transpose(Vec2);

M = [n(V0,V0), n(V0,V2), n(V0,V4); ...
    n(V2,V0), n(V2,V2), n(V2,V4); ...
    n(V4,V0), n(V4,V2), n(V4,V4)];
V = [n(V0,y); n(V2,y); n(V4,y)];

tv = M\V;
a0 = tv(1);
a2 = tv(2);
a4 = tv(3);

end


function [alpha, beta] = fn_fit(v2, v4, v)

m = [v2*transpose(v2), v2*transpose(v4); ...
    v4*transpose(v2), v4*transpose(v4)];

RHS = [v2*transpose(v); v4*transpose(v)];

tv = m\RHS;
alpha = tv(1);
beta = tv(2);

end


function res = fn_norm(vx)

res = mean(abs(reshape(vx, [1, numel(vx)])));

end


function vx_nt = fn_nt(vx)

vx_nt = zeros(1,length(vx));
vx_len = vx(2:end) - vx(1:end-1);
vx_nt(1:end-1) = 1/2*vx_len;
vx_nt(2:end) = vx_nt(2:end) + 1/2*vx_len;

end
